/*
 * Routines for calculating the force (F) and its derivate (FDOT).
 * We have several versions, including:
 * 
 *  - Pure serial computation
 *  - CUDA implementation
 * 
 * the last one is in fact an interface to the real CUDA code, which
 * is much more complex and requires an special compiler provided by 
 * NVIDIA (nvcc) 
 *
 * TODO: add automatic detection of CUDA capabilities
 *
 */
 
   #include <stdio.h>
   #include <stdlib.h>
   #include <omp.h>
   #include <math.h>
   
   #include <common.h>
 
// Pure serial computation
  void FFDOT( long i, double *F_i, double *FDOT_i, SIM_DATA *s )
  {

    int j,k;
      
    double DR[3];
    double DV[3];
    double RIJ2,DR3I,DRDV;
    
    for( k=0; k<3; k++ )
    {
      F_i[k] = 0.0;
      FDOT_i[k] = 0.0;
    }
    
    for( j=0; j< s->N; j++ )
    if(  j != i )
    {
      for( k=0; k<3; k++ )
      {
        DR[k] =    s->X[4*j+k] - s->X[4*i+k];
        DV[k] = s->XDOT[4*j+k] - s->XDOT[4*i+k];
      }
      
      RIJ2 = DR[0]*DR[0] + DR[1]*DR[1] + DR[2]*DR[2] + s->eps2;
      DR3I = s->X[4*j+3]/( RIJ2 * sqrt(RIJ2) );
      
      DRDV = 3.0*( DR[0]*DV[0] + DR[1]*DV[1] + DR[2]*DV[2] )/RIJ2;
      
      for( k=0; k<3; k++ )
      {
        F_i[k]    +=   DR[k]*DR3I;
        FDOT_i[k] += ( DV[k] - DR[k]*DRDV )*DR3I;
      }  
    }
  
  }

  double energy_check( SIM_DATA *s ) 
  {
    
    long i,j,tid;
    
    double K = 0.0,U = 0.0;
    
    double U_t[256];
    for( i=0; i<256; i++ )
		U_t[i] = 0.0;
		
    double *X = s->X;
    double *XDOT = s->XDOT;
    long Np = s->N;
  
    for( i=0; i<Np; i++ )
    {
      K += X[4*i+3]*(pow(XDOT[4*i+0],2)+pow(XDOT[4*i+1],2)+pow(XDOT[4*i+2],2));
      
      #pragma omp parallel for private(j,U_t,tid)
      for( j=i+1; j<Np; j++ )
	  {
        tid = omp_get_thread_num();
        
        
        U_t[tid] += X[4*i+3]*X[4*j+3]/sqrt( pow(X[4*j+0]-X[4*i+0],2) +
                                            pow(X[4*j+1]-X[4*i+1],2) + 
                                            pow(X[4*j+2]-X[4*i+2],2) + s->eps2 ); 
	  }
	  
	  for( j=0; j<omp_get_num_threads(); j++ )
	    U += U_t[j];
    }
    
    return 0.5*K - U;
  }
  
